# -*- coding: utf-8 -*-
"""
Created on Mon Nov 18 13:00:38 2013

@author: wudong

File location :
'D:/MyCodeLib/MyAbaqusLib/'

usage :
import numpy as np
import mdbHelper as mh

Reload the module :
import mdbHelper as mh
reload(mh)

"""
from abaqusConstants import *
import regionToolset
import numpy as np
import abaHelper as ah
#
# ============================================================================
#
def inRange(x, xRange, margin):
	(x1, x2) = xRange
	if x1>x2:
		xtemp = x2
		x2 = x1
		x1 = xtemp
	assert x1<=x2
	return ((x1-margin)<x) and (x<(x2+margin))
#
# ============================================================================
#
def getPartNodeLabelFromCoord(part, nodeCoord, tolerance=0.1):
	coord1 = nodeCoord
	if len(coord1)<3:
		coord1 = (nodeCoord[0], nodeCoord[1], 0)
	#
	allNodes = part.nodes
	nodeLabel = -1
	(x1,y1,z1) = coord1
	for node in allNodes:
		(x2,y2,z2) = node.coordinates
		dist = pow(x1-x2,2)+pow(y1-y2,2)+pow(z1-z2,2)
		if dist<tolerance:
			nodeLabel = node.label
			break
	assert (nodeLabel>=0)
	return nodeLabel
#
# ============================================================================
#
def getPartNodeLabelsFromCoords(part, nodeCoords, tolerance=0.1):
	nodeLabels = []
	for nodeCoord in nodeCoords:
		nodeLabel = getPartNodeLabelFromCoord(part, nodeCoord, tolerance)
		nodeLabels.append(nodeLabel)
	return nodeLabels
#
# ============================================================================
#
def createPartNodeSetFromLabels(part, nodeSetName, nodeLabels):
	nodeSet = part.SetFromNodeLabels(name=nodeSetName, nodeLabels=tuple(nodeLabels))
	return nodeSet
#
# ============================================================================
#
def getPartNodeLabelsFromNodeSet(part, nodeSetName):
	allNodes = part.sets[nodeSetName].nodes
	nodeLabels = []
	for node in allNodes:
		nodeLabels.append(node.label)
	return nodeLabels
#
# ============================================================================
#
def createPartNodeSetFromCoords(part, nodeSetName, nodeCoords, tolerance=0.1):
	nodeLabels = getPartNodeLabelsFromCoords(part, nodeCoords, tolerance)
	nodeSet = createPartNodeSetFromLabels(part, nodeSetName, nodeLabels)
	return nodeSet
#
# ============================================================================
#
def createPartNodeSetFromRectRange(part, nodeSetName, xRange, yRange, zRange=(0,0), margin = 0.1):
	allNodes = part.nodes
	nodeLabels = []
	for node in allNodes:
		(x,y,z) = node.coordinates
		if (inRange(x,xRange,margin) and inRange(y,yRange,margin) and inRange(z,zRange,margin)):
			nodeLabels.append(node.label)
	nodeSet = createPartNodeSetFromLabels(part, nodeSetName, nodeLabels)
	return nodeSet
#
# ============================================================================
#
def createPartElemSetFromLabels(part, elemSetName, elemLabels):
	elemSet = part.SetFromElementLabels(name=elemSetName, elementLabels=tuple(elemLabels))
	return elemSet
#
# ============================================================================
#
def createPartElemSetFromRectRange(part, elemSetName, xRange, yRange, zRange=(0,0), margin=0.1):
	allElems = part.elements
	elemLabels = []
	for elem in allElems:
		elemNodes = elem.getNodes()
		allNodesInRange = True
		for node in elemNodes:
			(x,y,z) = node.coordinates
			if not (inRange(x,xRange,margin) and inRange(y,yRange,margin) and inRange(z,zRange,margin)):
				allNodesInRange = False
				break
		if allNodesInRange:
			elemLabels.append(elem.label)
	elemSet = createPartElemSetFromLabels(part, elemSetName, elemLabels)
	return elemSet
#
# ============================================================================
#
#
# ============================================================================
#
def createPartNodeSpringDashpot(part, nodeLabel, dof, springStiffness=0, dashpotCoefficient=0, springDashpotName=None):
	assert (not ((springStiffness<=0) and (dashpotCoefficient<=0)))
	springBehavior = ON
	if springStiffness<=0:
		springBehavior = OFF
	#
	dashpotBehavior = ON
	if dashpotCoefficient<=0:
		dashpotBehavior = OFF
	#
	sdName = ''
	if springDashpotName==None:
		if springBehavior==ON:
			sdName = sdName+'Spring'
		if dashpotBehavior==ON:
			sdName = sdName+'Dashpot'
	sdName = '%s_%d_%d' % (sdName,nodeLabel,dof)
	#
	nodeRegion = part.nodes.sequenceFromLabels((nodeLabel,))
	region=regionToolset.Region(nodes=nodeRegion)
	part.engineeringFeatures.SpringDashpotToGround(
		name=sdName, region=region, orientation=None, dof=dof,
		springBehavior=springBehavior, springStiffness=springStiffness,
		dashpotBehavior=dashpotBehavior, dashpotCoefficient=dashpotCoefficient)
#
# ============================================================================
#
def createPartNodeSetSpringDashpot(part, nodeSetName, dof, springStiffness=0, dashpotCoefficient=0):
	assert (nodeSetName in part.sets.keys())
	allNodes = part.sets[nodeSetName].nodes
	for node in allNodes:
		nodeLabel = node.label
		createPartNodeSpringDashpot(part, nodeLabel, dof, springStiffness, dashpotCoefficient)
#
# ============================================================================
#
def getPartElemLabelsFromElemSet(part, elemSetName):
	allElems = part.sets[elemSetName].elements
	elemLabels = []
	for elem in allElems:
		elemLabels.append(elem.label)
	return elemLabels
#
# ============================================================================
#
def applyModelDynamicConcentratedForce(model, createStepName, nodeLabel, dof, ampData, bcName=None):
	assert (1<=dof) and (dof<=3)
	instance = model.rootAssembly.instances.values()[0]
	nodeRegion = instance.nodes.sequenceFromLabels((nodeLabel,))
	region=regionToolset.Region(nodes=nodeRegion)
	#
	_bcName = bcName
	if _bcName==None:
		_bcName = 'BC_CF_%s_%d_%d' % (createStepName, nodeLabel, dof)
	#
	ampName = 'Amp_%s' % (_bcName, )
	model.TabularAmplitude(name=ampName, timeSpan=STEP,
		smooth=SOLVER_DEFAULT, data=ah.arrayToTable(ampData))
	#
	cf1 = 0
	cf2 = 0
	cf3 = 0
	if dof==1:
		cf1 = 1.0
	elif dof==2:
		cf2 = 1.0
	else:
		cf3 = 1.0
	model.ConcentratedForce(name=_bcName, createStepName=createStepName,
		region=region, cf1=cf1, cf2=cf2, cf3=cf3, amplitude=ampName,
		distributionType=UNIFORM, field='', localCsys=None)
#
# ============================================================================
#
def applyModelDynamicDisplacement(model, createStepName, nodeLabel, dof, ampData, bcName=None):
	assert (1<=dof) and (dof<=3)
	instance = model.rootAssembly.instances.values()[0]
	nodeRegion = instance.nodes.sequenceFromLabels((nodeLabel,))
	region=regionToolset.Region(nodes=nodeRegion)
	#
	_bcName = bcName
	if _bcName==None:
		_bcName = 'BC_%s_%d_%d' % (createStepName, nodeLabel, dof)
	#
	ampName = 'Amp_%s' % (_bcName, )
	model.TabularAmplitude(name=ampName, timeSpan=STEP,
		smooth=SOLVER_DEFAULT, data=ah.arrayToTable(ampData))
	#
	u1 = UNSET
	u2 = UNSET
	u3 = UNSET
	if dof==1:
		u1 = 1.0
	elif dof==2:
		u2 = 1.0
	else:
		u3 = 1.0
	model.DisplacementBC(name=_bcName, createStepName=createStepName,
		region=region, u1=u1, u2=u2, ur3=u3, amplitude=ampName,
		distributionType=UNIFORM, fieldName='', localCsys=None)
#
# ============================================================================
#
def createModelNodeSetHistoryOutput(model, createStepName, nodeSetName, outputVarNames, numIntervals=500, histOutputName=None):
	instance = model.rootAssembly.instances.values()[0]
	assert (nodeSetName in instance.sets.keys())
	name = histOutputName
	if name==None:
		name = 'Hist_%s_%s' % (nodeSetName, '_'.join(outputVarNames))
	regionDef=instance.sets[nodeSetName]
	model.HistoryOutputRequest(name=name, createStepName=createStepName,
		variables=outputVarNames, region=regionDef,
		numIntervals=numIntervals,
		sectionPoints=DEFAULT, rebar=EXCLUDE)
#
# ============================================================================
#
def createModelRegionFieldOutput(model, createStepName, outputVarNames, numIntervals=100, region=MODEL, fieldOutputName=None):
	name = fieldOutputName
	if name==None:
		name = 'Field_region_%s' % ('_'.join(outputVarNames),)
	model.FieldOutputRequest(name=name, createStepName=createStepName,
		variables=outputVarNames,
		region=region, numIntervals=numIntervals)
#
# ============================================================================
# ============================================================================
#
from abaqus import *
from abaqusConstants import *
from caeModules import *
#
class AbaqusMdb:
	def __init__(self, modelName, inpFileName):
		self.modelName = modelName
		self.mdb = Mdb()
		self.model = self.mdb.ModelFromInputFile(name=self.modelName, 
										   inputFileName=inpFileName)
		self.partName = 'PART_ANSMODEL'
		self.part = self.model.parts[self.partName]
		pass

	#
	# material
	#
	def createElasticMaterial(self, matName, density, elasticModulus, poissonRatio):
		assert (0<density) and (density<10000)
		assert elasticModulus>0
		assert (0<poissonRatio) and (poissonRatio<0.5)
		mat = self.model.Material(name=matName)
		mat.Density(table=((density, ), ))
		mat.Elastic(table=((elasticModulus, poissonRatio), ))
		return mat

	def createElasticPlasticMaterial(self, matName, density, elasticModulus, poissonRatio, plasticData):
		assert (0<density) and (density<10000)
		assert elasticModulus>0
		assert (0<poissonRatio) and (poissonRatio<0.5)
		assert type(plasticData) is tuple
		mat = self.model.Material(name=matName)
		mat.Density(table=((density, ), ))
		mat.Elastic(table=((elasticModulus, poissonRatio), ))
		mat.Plastic(table=plasticData)
		return mat

	def createConcreteDamagedPlasticityMaterial(self, matName, density, elasticModulus, poissonRatio, 
				chData, cdData, tsData, tdData):
		mat = self.model.Material(name=matName)
		mat.Density(table=((density, ), ))
		mat.Elastic(table=((elasticModulus, poissonRatio), ))
		mat.ConcreteDamagedPlasticity(table=((15.0, 0.1, 1.16, 0.6667, 5e-4), ))
		cdp = mat.concreteDamagedPlasticity
		cdp.ConcreteCompressionHardening(table=chData)
		cdp.ConcreteCompressionDamage(table=cdData)
		#cdp.ConcreteTensionStiffening(table=tsData, type=DISPLACEMENT)
		#cdp.ConcreteTensionDamage(table=tdData, type=DISPLACEMENT)
		cdp.ConcreteTensionStiffening(table=tsData)
		cdp.ConcreteTensionDamage(table=tdData)
		return mat

	def assignMaterial(self, matName, esetName, thickness=None):
		assert esetName.upper() in self.part.sets.keys(), 'esetName=%s' % (esetName,)
		sectionName = 'Section_%s_%s' % (matName, esetName)
		self.model.HomogeneousSolidSection(name=sectionName, material=matName, thickness=thickness)
		region = self.part.sets[esetName.upper()]
		self.part.SectionAssignment(region=region, sectionName=sectionName, offset=0.0,
			offsetType=MIDDLE_SURFACE, offsetField='', thicknessAssignment=FROM_SECTION)
		pass

	def assignMaterialToSolidElemSet(self, matName, esetName, thickness=None):
		assert esetName.upper() in self.part.sets.keys(), 'esetName=%s' % (esetName,)
		sectionName = 'Section_%s_%s' % (matName, esetName)
		self.model.HomogeneousSolidSection(name=sectionName, material=matName, thickness=thickness)
		region = self.part.sets[esetName.upper()]
		self.part.SectionAssignment(region=region, sectionName=sectionName, offset=0.0,
			offsetType=MIDDLE_SURFACE, offsetField='', thicknessAssignment=FROM_SECTION)
		pass
	
	def assignMaterialToTrussElemSet(self, matName, esetName, sectionArea):
		assert esetName.upper() in self.part.sets.keys(), 'esetName=%s' % (esetName,)
		sectionName = 'Section_%s_%s' % (matName, esetName)
		self.model.TrussSection(name=sectionName, material=matName, area=sectionArea)
		region = self.part.sets[esetName.upper()]
		self.part.SectionAssignment(region=region, sectionName=sectionName, offset=0.0,
			offsetType=MIDDLE_SURFACE, offsetField='', thicknessAssignment=FROM_SECTION)
		pass

	def assignMaterialToRectBeamElemSet(self, matName, esetName, beamWidth, beamHeight, n1=(0.0, 0.0, -1.0)):
		assert esetName.upper() in self.part.sets.keys(), 'esetName=%s' % (esetName,)
		profileName = 'Profile_%s_%s' % (matName, esetName)
		self.model.RectangularProfile(name=profileName, a=beamWidth, b=beamHeight)
		sectionName = 'Section_%s_%s' % (matName, esetName)
		self.model.BeamSection(name=sectionName, integration=DURING_ANALYSIS, 
						 poissonRatio=0.0, profile=profileName, material=matName, 
						 temperatureVar=LINEAR, consistentMassMatrix=False)
		region = self.part.sets[esetName.upper()]
		self.part.SectionAssignment(region=region, sectionName=sectionName, offset=0.0, 
			offsetType=MIDDLE_SURFACE, offsetField='', thicknessAssignment=FROM_SECTION)
		self.part.assignBeamSectionOrientation(region=region, method=N1_COSINES, n1=n1)
		pass

	#
	# instance
	#
	def createInstance(self):
		self.instanceName = self.partName+'-1'
		self.instance = self.model.rootAssembly.Instance(name=self.instanceName, 
												   part=self.part, dependent=ON)
		pass

	#
	# create analysis step
	#
	def createStaticStep(self, stepName, timePeriod=1.0, previous='Initial', 
				initialInc=None, maxInc=None, minInc=None, maxNumInc=None, 
				nlgeom=OFF):
		assert timePeriod>0
		if maxInc==None:
			maxInc = timePeriod
		if initialInc==None:
			initialInc = maxInc
		if minInc==None:
			minInc = 1.0e-05
		if maxNumInc==None:
			maxNumInc = 100000
		self.model.StaticStep(name=stepName, previous=previous, 
				timePeriod=timePeriod, maxNumInc=maxNumInc, initialInc=initialInc, 
				minInc=minInc, maxInc=maxInc, nlgeom=nlgeom)
		pass

	def createImplicitDynamicsStep(self, stepName, timePeriod, previous='Initial', 
				initialInc=None, maxInc=None, minInc=None, maxNumInc=None, 
				nlgeom=OFF, matrixStorage=SOLVER_DEFAULT):
		assert timePeriod>0
		if maxInc==None:
			maxInc = timePeriod
		if initialInc==None:
			initialInc = maxInc
		if minInc==None:
			minInc = 1.0e-05
		if maxNumInc==None:
			maxNumInc = 100000
		self.model.ImplicitDynamicsStep(name=stepName, previous=previous, 
				timePeriod=timePeriod, maxNumInc=maxNumInc, initialInc=initialInc, 
				minInc=minInc, maxInc=maxInc, nlgeom=nlgeom, matrixStorage=matrixStorage)
		pass
		
	#
	# define boundary condition
	#
	def createGravityLoad(self, createStepName, comp1=UNSET, comp2=UNSET, comp3=UNSET):
		assert ((comp1!=UNSET and comp2==UNSET and comp3==UNSET) 
			or (comp1==UNSET and comp2!=UNSET and comp3==UNSET) 
			or (comp1==UNSET and comp2==UNSET and comp3!=UNSET))
		_bcName = 'BC_%s_Gravity' % (createStepName, )
		self.model.Gravity(name=_bcName, createStepName=createStepName, 
					 comp1=comp1, comp2=comp2, comp3=comp3,
					 distributionType=UNIFORM, field='')
		pass

	def importInitialStressFieldFromOdbFile(self, odbFileName, step=1, increment=1):
		_bcName = 'BC_PredefinedStressField_'
		self.model.Stress(name=_bcName, distributionType=FROM_FILE, fileName=odbFileName, 
						step=step, increment=increment)
		pass

	def createFixBC(self, createStepName, nsetName, fixU1=True, fixU2=True, fixU3=True, bcName=None):
		assert nsetName.upper() in self.instance.sets.keys()
		region = self.instance.sets[nsetName.upper()]
		_bcName = bcName
		if _bcName==None:
			_bcName = 'BC_Fix_%s_%s' % (createStepName, nsetName)
			if fixU1:
				_bcName = _bcName+'_U1'
			if fixU2:
				_bcName = _bcName+'_U2'
			if fixU3:
				_bcName = _bcName+'_U3'
		#
		u1Set = SET if fixU1 else UNSET
		u2Set = SET if fixU2 else UNSET
		u3Set = SET if fixU3 else UNSET
		self.model.DisplacementBC(name=_bcName, createStepName=createStepName, 
			region=region, u1=u1Set, u2=u2Set, ur3=u3Set, amplitude=UNSET, 
			distributionType=UNIFORM, fieldName='', localCsys=None)
		pass

	def createStaticDisplacement(self, createStepName, nodeName, dof, dispAmp, bcName=None):
		assert (1<=dof) and (dof<=3)
		assert (type(nodeName) is str) or (type(nodeName) is int)
		region = None
		if type(nodeName) is str:		# nodeName is nodeSetName
			assert nodeName.upper() in self.instance.sets.keys()
			region = self.instance.sets[nodeName.upper()]
		elif type(nodeName) is int:		# nodeName is nodeLabel
			nodeRegion = self.instance.nodes.sequenceFromLabels((nodeName,))
			region=regionToolset.Region(nodes=nodeRegion)
		assert region != None
		#
		_bcName = bcName
		if _bcName==None:
			_bcName = 'BC_Disp_%s_%s_%d' % (createStepName, str(nodeName), dof)
		
		#
		u1 = UNSET
		u2 = UNSET
		u3 = UNSET
		if dof==1:
			u1 = dispAmp
		elif dof==2:
			u2 = dispAmp
		else:
			u3 = dispAmp
		self.model.DisplacementBC(name=_bcName, createStepName=createStepName,
			region=region, u1=u1, u2=u2, ur3=u3,
			distributionType=UNIFORM, fieldName='', localCsys=None)
		pass

	def createDynamicConcentratedForce(self, createStepName, nodeName, dof, ampData, bcName=None):
		assert (1<=dof) and (dof<=3)
		assert (type(nodeName) is str) or (type(nodeName) is int)
		assert (type(ampData) is np.ndarray), 'type(ampData) is %s'%(str(type(ampData)))
		assert (len(ampData.shape)==2) and (ampData.shape[1]==2), 'ampData.shape = %s' % (str(ampData.shape))
		region = None
		if type(nodeName) is str:		# nodeName is nodeSetName
			assert nodeName.upper() in self.instance.sets.keys()
			region = self.instance.sets[nodeName.upper()]
		elif type(nodeName) is int:		# nodeName is nodeLabel
			nodeRegion = self.instance.nodes.sequenceFromLabels((nodeName,))
			region=regionToolset.Region(nodes=nodeRegion)
		assert region != None
		#
		_bcName = bcName
		if _bcName==None:
			_bcName = 'BC_CF_%s_%s_%d' % (createStepName, nodeName, dof)
		#
		ampName = 'Amp_%s' % (_bcName, )
		self.model.TabularAmplitude(name=ampName, timeSpan=STEP,
			smooth=SOLVER_DEFAULT, data=ah.arrayToTable(ampData))
		#
		cf1 = 0
		cf2 = 0
		cf3 = 0
		if dof==1:
			cf1 = 1.0
		elif dof==2:
			cf2 = 1.0
		else:
			cf3 = 1.0
		self.model.ConcentratedForce(name=_bcName, createStepName=createStepName,
			region=region, cf1=cf1, cf2=cf2, cf3=cf3, amplitude=ampName,
			distributionType=UNIFORM, field='', localCsys=None)
		pass

	def createNodeSetDynamicConcentratedForce(self, createStepName, nodeData, dof, ampData):
		assert (1<=dof) and (dof<=3)
		assert nodeData.shape[1]==5		# nodeData is a numpy array exported by '~exni' command from ANSYS
		nodeCount = nodeData.shape[0]
		assert nodeCount>0
		_ampData = np.zeros(shape=ampData.shape)
		_ampData[:, 0] = ampData[:, 0]	# t
		for i in range(nodeCount):
			nodeLabel = int(nodeData[i, 0])
			nodeArea = nodeData[i, 4]
			_ampData[:, 1] = ampData[:, 1]*nodeArea
			self.createDynamicConcentratedForce(createStepName, nodeLabel, dof, _ampData)
		pass

	def createDynamicDisplacement(self, createStepName, nodeName, dof, ampData, bcName=None):
		assert (1<=dof) and (dof<=3)
		assert (type(nodeName) is str) or (type(nodeName) is int)
		region = None
		if type(nodeName) is str:		# nodeName is nodeSetName
			assert nodeName.upper() in self.instance.sets.keys()
			region = self.instance.sets[nodeName.upper()]
		elif type(nodeName) is int:		# nodeName is nodeLabel
			nodeRegion = self.instance.nodes.sequenceFromLabels((nodeName,))
			region=regionToolset.Region(nodes=nodeRegion)
		assert region != None
		#
		_bcName = bcName
		if _bcName==None:
			_bcName = 'BC_Disp_%s_%s_%d' % (createStepName, str(nodeName), dof)
		#
		ampName = 'Amp_%s' % (_bcName, )
		self.model.TabularAmplitude(name=ampName, timeSpan=STEP,
			smooth=SOLVER_DEFAULT, data=ah.arrayToTable(ampData))
		#
		u1 = UNSET
		u2 = UNSET
		u3 = UNSET
		if dof==1:
			u1 = 1.0
		elif dof==2:
			u2 = 1.0
		else:
			u3 = 1.0
		self.model.DisplacementBC(name=_bcName, createStepName=createStepName,
			region=region, u1=u1, u2=u2, ur3=u3, amplitude=ampName,
			distributionType=UNIFORM, fieldName='', localCsys=None)
		pass

	def createNodeSpringDashpot(self, nodeLabel, dof, springStiffness=0, dashpotCoefficient=0, springDashpotName=None):
		assert ((springStiffness>0) or (dashpotCoefficient>0))
		springBehavior = ON if springStiffness>0 else OFF
		dashpotBehavior = ON if dashpotCoefficient>0 else OFF
		#
		sdName = ''
		if springDashpotName==None:
			if springBehavior==ON:
				sdName = sdName+'Spring'
			if dashpotBehavior==ON:
				sdName = sdName+'Dashpot'
		sdName = '%s_%d_%d' % (sdName,nodeLabel,dof)
		#
		nodeRegion = self.part.nodes.sequenceFromLabels((nodeLabel,))
		region=regionToolset.Region(nodes=nodeRegion)
		self.part.engineeringFeatures.SpringDashpotToGround(
			name=sdName, region=region, orientation=None, dof=dof,
			springBehavior=springBehavior, springStiffness=springStiffness,
			dashpotBehavior=dashpotBehavior, dashpotCoefficient=dashpotCoefficient)
		pass

	def createNodeSetSpringDashpot(self, nodeData, dof, Kb=0, Cb=0):
		assert Kb>0 or Cb>0
		assert nodeData.shape[1]==5		# nodeData is a numpy array exported by '~exni' command from ANSYS
		nodeCount = nodeData.shape[0]
		assert nodeCount>0
		for i in range(nodeCount):
			nodeLabel = int(nodeData[i, 0])
			nodeArea = nodeData[i, 4]
			springStiffness = Kb*nodeArea
			dashpotCoefficient = Cb*nodeArea
			self.createNodeSpringDashpot(nodeLabel, dof, springStiffness=springStiffness, 
								dashpotCoefficient=dashpotCoefficient)
		pass

	def applyDRMNodesPeb(self, createStepName, nodeLabels, timePointData, Peb):
		nodeCount = len(nodeLabels)
		totalDofCount = Peb.shape[0]
		assert (totalDofCount % nodeCount)==0
		dofCount = totalDofCount/nodeCount
		timePointCount = timePointData.shape[0]
		assert timePointCount==Peb.shape[1]
		#
		ampData = np.zeros(shape=(timePointCount, 2))
		ampData[:,0] = timePointData.reshape((timePointCount,))
		for i in range(nodeCount):
			nodeLabel = nodeLabels[i]
			for j in range(dofCount):
				ampData[:,1] = Peb[i*dofCount+j,:].reshape((timePointCount,))
				dof = j+1
				self.createDynamicConcentratedForce(createStepName, nodeLabel, dof, ampData)
		pass

	#
	# define output request
	#
	def createHistoryOutput(self, histOutputName, createStepName, variables, setName, 
						 numIntervals=100):
		assert setName.upper() in self.instance.sets.keys()
		region=self.instance.sets[setName.upper()]
		self.model.HistoryOutputRequest(name=histOutputName, createStepName=createStepName, 
			variables=variables, region=region, numIntervals=numIntervals, 
			sectionPoints=DEFAULT, rebar=EXCLUDE)
		pass

	def createFieldOutput(self, fieldOutputName, createStepName, variables, 
					   esetName=None, nsetName=None, numIntervals=10):
		region = MODEL
		eset = None
		nset = None
		if esetName != None :
			assert esetName.upper() in self.instance.sets.keys()
			eset = self.instance.sets[esetName.upper()].elements
		if nsetName != None :
			assert nsetName.upper() in self.instance.sets.keys()
			nset = self.instance.sets[nsetName.upper()].nodes
		if eset!=None or nset!=None :
			region = regionToolset.Region(elements=eset, nodes=nset)
		#
		self.model.FieldOutputRequest(name=fieldOutputName, createStepName=createStepName, 
			variables=variables, region=region, numIntervals=numIntervals)
		pass

	#
	# save model
	#
	def saveModel(self, caeFileName):
		self.mdb.saveAs(pathName=caeFileName)
		pass

	#
	# create/submit job
	#
	def createJob(self, jobName, cpuCount = 1):
		assert cpuCount>0
		self.job = self.mdb.Job(name=jobName, model=self.modelName, 
						  numCpus=cpuCount, numDomains=cpuCount)
		pass

	def submitJob(self):
		self.job.writeInput(consistencyChecking=OFF)
		self.job.submit(consistencyChecking=OFF)
		pass
	
	def waitForJobCompletion(self):
		self.job.waitForCompletion()
		pass

	#
	# insert/append additional input commands
	#
	def getModelBlockPosition(self, blockPrefix):
		if blockPrefix == '':
			return len(self.model.keywordBlock.sieBlocks)-1
		pos = 0
		for block in self.model.keywordBlock.sieBlocks:
			if block[0:len(blockPrefix)].lower()==blockPrefix.lower():
				return pos
			pos=pos+1
		return -1
	
	def insertMkMatrixOutputCommand(self, infElemSetName):
		self.model.keywordBlock.synchVersions()
		blockPrefix = '*End Step'
		insertPos = self.getModelBlockPosition(blockPrefix)-1
		assert insertPos>=0
		mkDataOuputStr = (
"""**
** Output Mass and stiffness matrix
**
*ELEMENT MATRIX OUTPUT, ELSET=%s.%s, MASS=YES, STIFFNESS=YES, FREQUENCY=1, OUTPUT FILE=USER DEFINED, FILE NAME=MK_data.txt
**""" % (self.instanceName, infElemSetName, ))
		self.model.keywordBlock.insert(insertPos, mkDataOuputStr)
		pass
	
	#
	# get/create node set
	#
	def getPartNodeLabelFromCoord(self, nodeCoord, tolerance=0.1):
		coord1 = nodeCoord
		if len(coord1)<3:
			coord1 = (nodeCoord[0], nodeCoord[1], 0)
		#
		allNodes = self.part.nodes
		nodeLabel = -1
		(x1,y1,z1) = coord1
		for node in allNodes:
			(x2,y2,z2) = node.coordinates
			dist = sqrt(pow(x1-x2,2)+pow(y1-y2,2)+pow(z1-z2,2))
			if dist<tolerance:
				nodeLabel = node.label
				break
		assert (nodeLabel>=0)
		return nodeLabel

	def getPartNodeLabelsFromCoords(self, nodeCoords, tolerance=0.1):
		nodeLabels = []
		for nodeCoord in nodeCoords:
			nodeLabel = self.getPartNodeLabelFromCoord(nodeCoord, tolerance)
			nodeLabels.append(nodeLabel)
		return nodeLabels

	def getPartNodeLabelsFromNodeSet(self, nodeSetName):
		assert nodeSetName.upper() in self.part.sets.keys()
		allNodes = self.part.sets[nodeSetName.upper()].nodes
		nodeLabels = [node.label for node in allNodes]
		return nodeLabels

	def createPartNodeSetFromLabels(self, nodeSetName, nodeLabels):
		nodeSet = self.part.SetFromNodeLabels(name=nodeSetName, nodeLabels=tuple(nodeLabels))
		return nodeSet

	def createPartNodeSetFromCoords(self, nodeSetName, nodeCoords, tolerance=0.1):
		nodeLabels = self.getPartNodeLabelsFromCoords(nodeCoords, tolerance)
		nodeSet = self.createPartNodeSetFromLabels(nodeSetName, nodeLabels)
		return nodeSet

	def getPartNodeCoordinateFromNodeLabel(self, nodeLabel):
		allNodes = self.part.nodes
		nodeCoord = None
		for node in allNodes:
			if node.label == nodeLabel:
				nodeCoord = node.coordinates
		assert nodeCoord!=None
		return nodeCoord

	#
	# create surface object from element set and node set
	#
	def createSurfaceFromElemSetAndNodeSet(self, surfaceName, elemSetName, nodeSetName):
		assert elemSetName.upper() in self.part.sets.keys()
		assert nodeSetName.upper() in self.part.sets.keys()
		elems = self.part.sets[elemSetName.upper()].elements
		nodes = self.part.sets[nodeSetName.upper()].nodes
		elemLabels = [elem.label for elem in elems]
		nodeLabels = [node.label for node in nodes]
		assert len([elem.label for elem in elems if len(elem.getNodes())!=4])<=0
		f1ElemLabels = [elem.label for elem in elems if ((elem.getNodes()[0].label in nodeLabels) and (elem.getNodes()[1].label in nodeLabels))]
		f2ElemLabels = [elem.label for elem in elems if ((elem.getNodes()[1].label in nodeLabels) and (elem.getNodes()[2].label in nodeLabels))]
		f3ElemLabels = [elem.label for elem in elems if ((elem.getNodes()[2].label in nodeLabels) and (elem.getNodes()[3].label in nodeLabels))]
		f4ElemLabels = [elem.label for elem in elems if ((elem.getNodes()[3].label in nodeLabels) and (elem.getNodes()[0].label in nodeLabels))]
		f1Elems = elems.sequenceFromLabels(f1ElemLabels)
		f2Elems = elems.sequenceFromLabels(f2ElemLabels)
		f3Elems = elems.sequenceFromLabels(f3ElemLabels)
		f4Elems = elems.sequenceFromLabels(f4ElemLabels)
		self.part.Surface(face1Elements=f1Elems, face2Elements=f2Elems, face3Elements=f3Elems, face4Elements=f4Elems, name=surfaceName)
		pass

	def createSurfaceToSurfaceContact(self, createStepName, masterSurfaceName, slaveSurfaceName, frictionFactor=0.1):
		assert masterSurfaceName in self.instance.surfaces.keys()
		assert slaveSurfaceName in self.instance.surfaces.keys()
		contactPropName = 'ContactProp_%s_%s' % (masterSurfaceName, slaveSurfaceName)
		contactProp = self.model.ContactProperty(contactPropName)
		contactProp.TangentialBehavior(formulation=PENALTY, directionality=ISOTROPIC, slipRateDependency=OFF, 
			pressureDependency=OFF, temperatureDependency=OFF, dependencies=0, table=((frictionFactor, ), ), 
			shearStressLimit=None, maximumElasticSlip=FRACTION, fraction=0.005, elasticSlipStiffness=None)
		contactProp.NormalBehavior(pressureOverclosure=HARD, allowSeparation=ON, constraintEnforcementMethod=DEFAULT)
		#
		region1=self.instance.surfaces[masterSurfaceName]
		region2=self.instance.surfaces[slaveSurfaceName]
		contactName = 'Contact_%s_%s' % (masterSurfaceName, slaveSurfaceName)
		self.model.SurfaceToSurfaceContactStd(name=contactName, 
			createStepName=createStepName, master=region1, slave=region2, sliding=SMALL, 
			thickness=ON, interactionProperty=contactPropName, 
			adjustMethod=NONE, initialClearance=OMIT, datumAxis=None, clearanceRegion=None)
		pass

	def createTieConstrainForNodeSets(self, masterNodeSetName, slaveNodeSetName):
		assert masterNodeSetName.upper() in self.instance.sets.keys()
		assert slaveNodeSetName.upper() in self.instance.sets.keys()
		region1=self.instance.sets[masterNodeSetName.upper()]
		region2=self.instance.sets[slaveNodeSetName.upper()]
		constrainName = 'Constrain_%s_%s' % (masterNodeSetName, slaveNodeSetName)
		self.model.Tie(name=constrainName, master=region1, slave=region2, positionToleranceMethod=COMPUTED, 
				 adjust=ON, tieRotations=ON, thickness=ON)
		pass

	#
